ct = data4d('../data/lung/ct')
itv = data3d('../data/lung/itv')
step = float_param(0.001, 0.0005, 0.01, 'S')

# rpl params
colormap = data1d_rgba('jet.png')
table = data1d('density-to-rpl-table.png')
sigma_max = float_param(0.1, 0, 0.5, 'M')

# transfer function params
tf = data1d_rgba('bone.png')
tf_pos = float_param(0, -200, 200, 'P')
tf_width = float_param(100, 10, 400, 'W')
tf_opacity = float_param(0.8, 0, 1, 'O')

# lighting params
ambient = float_param(0.1, 0, 1, 'A')
lightpos = float3_param(-3, 3, 3)

cut = float_param(-2, -2, 2, 'C')

C = normalize(S - E)
dP = (-C) * step
ray_length = length(E - S)
    
# first pass: find distal edge of ITV
t = 0.0
distal = -1.0
P = S

while t < ray_length:
    q = linear_query_3d(itv, P)
    if q > 0: distal = t
    P += dP
    t += step

# second pass: compute RPL sigma
sigma_color = float3(1)

if distal > 0:
    avg_squares = 0.0
    avg = 0.0

    for frame in range(10):
        P = S
        t = 0.0
        rpl = 0.0

        while t < distal:
            density = linear_query_4d(ct, P, frame) * 32768
            rel_density = 0.21 + (1.78 - 0.21) * linear_query_1d(table,
                (density - (-766)) / (1173 - (-766)))
            rpl += rel_density * step
            P += dP
            t += step
        
        avg += rpl
        avg_squares += rpl * rpl
        
    avg /= 10
    avg_squares /= 10
    rel_sigma = sqrt(avg_squares - avg * avg) / avg
    sigma_color = linear_query_1d_rgba(colormap, rel_sigma / sigma_max).xyz
    
# third pass: volume render
t = 0.0
P = S
result = float3(0)
see_through = 1.0

while t < ray_length:
    density = cubic_query_4d_cut(ct, P, C, cut, 0) * 32768
    tf_query = (density - tf_pos + tf_width) / (2 * tf_width)

    if tf_query >= 0:
        rgba = linear_query_1d_rgba(tf, tf_query)

        color = rgba.xyz * sigma_color
        
        N = -normalize(cubic_gradient_4d_cut(ct, P, C, cut, 0))
        L = normalize(lightpos - P)
        lit = phong(L, N, C, color, 0.5, 50, ambient)

        result += see_through * rgba.w * tf_opacity * lit
        see_through *= 1 - rgba.w * tf_opacity
        if see_through < 0.01: break

    P += dP
    t += step

return result
